######################################################################
### R file to examine out-of-sample prediction accuracy across weights matrices for  "All Economics Is Local"
###
### Created: 3-4-16
### Modified:
###
### NOTE: produces the bootstrapped percentiles in Table 3
######################################################################


## load packages and read data ##

library(foreign)
library(MASS)
library(mnormt)
library(nonnest2)

## Set up working directory ##
#setwd("")
	
data <- read.dta('All Economics Is Local.dta')

contNat <- c()
crossCont <- c()
cosponsorCross <- c()
econCosponsor <- c()
mediaEcon <- c()

set.seed(0624)

## estimate all the models then calculate predicted values ##

for(count in 1:1000){
	draw <- rank(runif(data$year))[1:(length(data$year)/5)]

	national <- polr(retnat ~ gsp_pc2_ch + gdppc_growth + ch_unem + inflation + approve + age +
		male + nonwhite + union + college + married + unemployed + own_home + ppid + noppid,
		data = data[-draw, ], method = 'logistic', Hess = TRUE, na.action = na.exclude)
	natVal <- predict(national, type = 'probs', na.action = na.exclude, newdata = data[draw, ])
	natPreds <- natVal[, 1] + (natVal[, 2] * 2) + (natVal[, 3] * 3)
	natDiff <- natPreds - as.numeric(data$retnat[draw])
	natDiff <- na.omit(natDiff)

	cont <- polr(retnat ~ gsp_pc2_ch + sp_W2_gsp_pc2_ch + ch_unem + inflation + approve + age +
		male + nonwhite + union + college + married + unemployed + own_home + ppid + noppid,
		data = data[-draw, ], method = 'logistic', Hess = TRUE, na.action = na.exclude)
	contVal <- predict(cont, type = 'probs', na.action = na.exclude, newdata = data[draw, ])
	contPreds <- contVal[, 1] + (contVal[, 2] * 2) + (contVal[, 3] * 3)
	contDiff <- contPreds - as.numeric(data$retnat[draw])
	contDiff <- na.omit(contDiff)

	cosponsor <- polr(retnat ~ gsp_pc2_ch + sp_W7_gsp_pc2_ch + ch_unem + inflation + approve + age +
		male + nonwhite + union + college + married + unemployed + own_home + ppid + noppid,
		data = data[-draw, ], method = 'logistic', Hess = TRUE, na.action = na.exclude)
	cosponsorVal <- predict(cosponsor, type = 'probs', na.action = na.exclude, newdata = data[draw, ])
	cosponsorPreds <- cosponsorVal[, 1] + (cosponsorVal[, 2] * 2) + (cosponsorVal[, 3] * 3)
	cosponsorDiff <- cosponsorPreds - as.numeric(data$retnat[draw])
	cosponsorDiff <- na.omit(cosponsorDiff)

	media <- polr(retnat ~ gsp_pc2_ch + sp_W10U_gsp_pc2_ch + ch_unem + inflation + approve + age +
		male + nonwhite + union + college + married + unemployed + own_home + ppid + noppid,
		data = data[-draw, ], method = 'logistic', Hess = TRUE, na.action = na.exclude)
	mediaVal <- predict(media, type = 'probs', na.action = na.exclude, newdata = data[draw, ])
	mediaPreds <- mediaVal[, 1] + (mediaVal[, 2] * 2) + (mediaVal[, 3] * 3)
	mediaDiff <- mediaPreds - as.numeric(data$retnat[draw])
	mediaDiff <- na.omit(mediaDiff)

	cross <- polr(retnat ~ gsp_pc2_ch + sp_W12_gsp_pc2_ch + ch_unem + inflation + approve + age +
		male + nonwhite + union + college + married + unemployed + own_home + ppid + noppid,
		data = data[-draw, ], method = 'logistic', Hess = TRUE, na.action = na.exclude)
	crossVal <- predict(cross, type = 'probs', na.action = na.exclude, newdata = data[draw, ])
	crossPreds <- crossVal[, 1] + (crossVal[, 2] * 2) + (crossVal[, 3] * 3)
	crossDiff <- crossPreds - as.numeric(data$retnat[draw])
	crossDiff <- na.omit(crossDiff)

	econ <- polr(retnat ~ gsp_pc2_ch + sp_W14_gsp_pc2_ch + ch_unem + inflation + approve + age +
		male + nonwhite + union + college + married + unemployed + own_home + ppid + noppid,
		data = data[-draw, ], method = 'logistic', Hess = TRUE, na.action = na.exclude)
	econVal <- predict(econ, type = 'probs', na.action = na.exclude, newdata = data[draw, ])
	econPreds <- econVal[, 1] + (econVal[, 2] * 2) + (econVal[, 3] * 3)
	econDiff <- econPreds - as.numeric(data$retnat[draw])
	econDiff <- na.omit(econDiff)

	contNat[count] <- 1 - pbinom(length(contDiff[abs(contDiff) < abs(natDiff)]), size = length(crossDiff), prob = 0.5)
	crossCont[count] <- 1 - pbinom(length(crossDiff[abs(crossDiff) < abs(contDiff)]), size = length(natDiff), prob = 0.5)
	cosponsorCross[count] <- 1 - pbinom(length(cosponsorDiff[abs(cosponsorDiff) < abs(crossDiff)]), size = length(cosponsorDiff), prob = 0.5)
	econCosponsor[count] <- 1 - pbinom(length(econDiff[abs(econDiff) < abs(cosponsorDiff)]), size = length(econDiff), prob = 0.5)
	mediaEcon[count] <- 1 - pbinom(length(mediaDiff[abs(mediaDiff) < abs(econDiff)]), size = length(mediaDiff), prob = 0.5)

	cat('run', count, 'of 1,000 complete...\n')
}

hold <- data.frame(contNat, crossCont, cosponsorCross, econCosponsor, mediaEcon)
write.table(hold, file = 'OOSPresults.txt', sep = '|')


#############################################################
## Display the 2.5, 50, and 97.5th percentiles

data <- read.table('OOSPresults.txt', header = TRUE, sep = '|')
summary(data)
	
meds <- c(median(data$contNat),
	median(data$crossCont), 
	median(data$cosponsorCross),
	median(data$econCosponsor),
	median(data$mediaEcon))

upper <- c(quantile(data$contNat, 0.975),
	quantile(data$crossCont, 0.975), 
	quantile(data$cosponsorCross, 0.975),
	quantile(data$econCosponsor, 0.975),
	quantile(data$mediaEcon, 0.975))

lower <- c(quantile(data$contNat, 0.025),
	quantile(data$crossCont, 0.025), 
	quantile(data$cosponsorCross, 0.025),
	quantile(data$econCosponsor, 0.025),
	quantile(data$mediaEcon, 0.025))

cbind(round(lower, 3), round(meds, 3), round(upper, 3))
	